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Abstract — We present a theoretical study on the 
dynamical properties of three-dimensional arrays of 
Josephson junctions. Our results indicate that such 
superconducting networks represent highly sensitive 
3D-SQUIDs having some major advantages in com- 
parison with conventional planar SQUIDs. The volt- 
age response function of 3D-SQUIDs is directly re- 
lated to the vector-character of external electromag- 
netic fields. The theory developed here allows the 
three-dimensional reconstruction of a detected exter- 
nal field including phase information about the field 
variables. Applications include the design of novel 
magnetometers, gradiometers and particle detectors. 



I. Introduction 

Oscillators based on arrays of Josephson junctions op- 
crating in the voltage mode have shown to be promis- 
ing radiation sources in the mm- and sub-mm wavelength 
range. One of the most remarkable features of such sys- 
tems is that the microwave power output increases with 
the number of array junctions while at the same time the 
linewidth of the generated high-frequency voltage oscilla- 
tion decreases with this number [6]. The latter phenom- 
ena is not only interesting for the application of the ar- 
rays as microwave sources but also for their application as 
highly sensitive detectors. From this point of view, two- 
and three-dimensional arrays of Josephson junctions rep- 
resent multi-junction interferometers which can be used 
e.g., as magnetometers, gradiometers or particle detec- 
tors. 

The capabilities of 2D- and 3D- Josephson arrays can 
be used successfully if all array junctions operate in a 
coherent array mode, i.e., if they oscillate strongly phase- 
locked [3,6]. The problem of phase-locking in 2D-arrays 
is presently subject of intensive research and great exper- 
imental and theoretical progress has been achieved using 
present-day integrated-circuit technology [1-3] . Since this 
technology is essentially two-dimensional there are up to 
now no corresponding experimental and theoretical stud- 
ies on 3D-arrays. However, 3D-arrays of Josephson junc- 
tions may be promising candidates for highly sensitive de- 
tectors which allow as a new quality the three-dimensional 
reconstruction of the detected electromagnetic field in- 
cluding phase information about the field variables. In 
view of recent developments in molecular beam epitaxy 
the fabrication of 3D-arrays with appropriate junction pa- 
rameters and array homogeneity seems to be possible in 
the nearest future. 
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Fig. 1. Schematic drawing of a 3D-network of Josephson junctions. 
The network is built by superconducting islands (large cubes) con- 
nected by Josephson junctions (smaller dashed cubes). 

In this paper we present theoretical results on the dy- 
namical properties of a 3D-system of coupled Joseph- 
son junctions. A schematic drawing of the 3D-network 
is shown in Fig.[j]. It consists of superconducting is- 
lands (large cubes) connected by Josephson junctions 
(smaller dashed cubes) . This system is the simplest three- 
dimensional configuration that guarantees coherent oper- 
ation, i.e., single frequency operation, if it is DC biased 
parallel to one of the networks' axes (c.f. Fig.||). If in 
the voltage state of the array a magnetic field is applied 
to the network, the six intrinsically coupled four-junction 
DC SQUIDs lying on the six faces of the cube are gen- 
erating field dependent macroscopic current and voltage 
distributions which can be easily measured. From the AC 
Josephson effect and the flux quantization condition the 
current and voltage distributions depend directly but non- 
linearly on the strength and the orientation of the external 
field. 



II. Network equations 

To get quantitative results about the networks' response 
function we will consider a network built by junctions 
which can be described by the RCSJ-model, e.g., exter- 
nally shunted Nb — A10 X ~ Nb junctions. However, if other 
SIS-type junctions are used the results remain qualita- 
tively valid. Fig|5 shows the equivalent circuit for the 3D- 
network which is~DC biased along the x-direction. It con- 
sists of Josephson junctions (crosses) which are connected 
by superconducting wires. The bias current Ib is fed into 
and extracted from the network through ohmic resistors 
(which must not necessarily be equal) . This method of bi- 
asing preserves well defined boundary conditions for the 
networks dynamics and removes some of the ambiguities 



1 




Fig. 2. Equivalent circuit of the 3D-network. The Josephson junc- 
tions (crosses) are connected by superconducting wires and the bias 
current Ig is fed into and extracted from the system through ohmic 
resistors (black rectangles). 

due to the symmetry of the system. Since the bias cur- 
rent is flowing parallel to the x-axis there are two groups 
of junctions with different dynamical behavior. The junc- 
tions lying parallel to the bias current (active junctions) 
switch into their voltage state for Ib > Ic,A(fext)> w here 
-^c,A(/ext) ^ s the critical current of the array depending 
on the externally applied flux / ex t = (fx, fy, /»)■ The 
junctions lying perpendicular to the bias current (passive 
junctions) do in general not switch into the voltage state 
but show librations (semirotations) of their Josephson- 
phase differences. The amplitude of these librations is 
directly related to the strength and orientation of the ex- 
ternal field. If the values of the ohmic resistors at the 
current in- and output do not differ too much, the pas- 
sive junctions arc in the zero-voltage state and the only 
frequency present in the system is the driving frequency 
defined by the active network junctions. This prevents 
chaotic or aperiodic dynamics of the junction phases. 

In the temporal gauge the variables characterizing the 
dynamics of the junction network are the gauge invariant 
phase differences [4,5]. If we denote the vector of the 12 
network variables by x(t) an d use the usual reduced units 
[1] the system of coupled network equations can be written 
symbolically as 

Px + X + sinx = A [L- l x + Ffexk] C 1 ) 
+Ej r{ 

+G(r B )I B 

+D(r B ) X , 

where A = &o/(2ir /j,ai c a) is the magnetic penetration 
depth, L the inductance matrix of the array, F the matrix 
of the distribution of the external flux / ex t over the differ- 
ent network meshes and E the matrix which describes the 
distribution of the currents j r j induced by the oscillating 
electric field vector of an incoming external wave. The el- 
ements of E and F depend on the orientation of the exter- 
nal electromagnetic field. G and D are the matrices which 
define the boundary conditions depending on the values 
of the eight input and output resistors r B — \TB,i} ■ The 
single junction parameters are the McCumber parameter 
(3, the critical current i c and the normal resistance r. <E>o 
is the magnetic flux quantum and a the lattice spacing 
measuring the distance between the centers of adjacent 
superconducting islands. The inductance matrix L in- 
cludes all inductances present in the array and depends 
on the geometry of the network and the inductances of 
the single junctions. We determine the coefficients of L 
by assuming an array geometry similar to that shown in 



Fig.|l| and by using junction inductances which are typical 
for externally shunted junctions. Changing the geometry 
of the network or the value of the junction inductances, 
however, influences our results only slightly. 

The network equation ([!]) clearly shows the dependence 
of the time evolution of the network phases x(i) on the 
applied quasi-static external magnetic field which can ap- 
proximately be expressed as / ex t = ^ext( x , Vi z ) a2 /^o, 
where -B ex t is the magnetic field vector of the external 
field. The sensitivity of the network depends on the lat- 
tice spacing a and reaches its maximum for vanishing mag- 
netic coupling within the network, i.e., for infinite mag- 
netic penetration depth A — > oo. For A — > 0, however, 
the network response vanishes totally. The sensitivity of 
the network dynamics with respect to distinct directions 
of the magnetic field, i.e., the anisotropy of the response 
function, can be manipulated by an appropriate choice of 
the input and output resistors. 

The device possesses in principle two modes of opera- 
tion. The DC mode with Ig < I c ^a and the AC mode 
with Ib > I C ,A- For subcritical bias currents Ib < I C ,A, a 
constant external field induces constant loop currents in 
the meshes of the network and the resulting current dis- 
tribution within the array is a superposition of these loop 
currents with the bias current. The critical array current 
-^c,A(/ext) ^ s a function of the strength of the external 
field and the orientation of this field. However, because of 
the many degrees of freedom of the system the response 
function is not unique. There exists in general for any 
given -B ex t a whole set of different equilibrium current 
distributions in the network on which the critical array 
current depends [5,7]. By a detailed theoretical analysis 
it is possible to formulate selection rules which describe 
the external field depending transitions between the dif- 
ferent distributions. Although there occur in the region of 
subcritical bias currents very interesting dynamical phe- 
nomena, the scenario becomes very complicated and will 
be presented elsewhere [7] . 

In the following we will restrict our treatment to the 
AC mode of the network. For Ib > 4i c the active ar- 
ray junctions show persistent voltage oscillations whose 
frequency is directly proportional to the voltage drop. 
If the constant voltage drop originating from the input 
and output resistors is subtracted, the averaged volt- 
age drop across the array is for Ib > I C ,A given by 
Vq = 1/4 < Ylx Xx(t) >, where the summation runs 
over the four active network junctions lying parallel to 
the x-axis and < . . . > denotes time averaging over one 
oscillation period. Due to flux quantization, the network 
currents induced by the external field are converted by the 
network dynamics into oscillating loop currents in each 
network mesh. The voltage oscillation induced by these 
oscillating loop currents interferes with the voltage oscil- 
lation of the active junctions, so that the averaged voltage 
drop V across the whole array becomes a function of the 
strength and the orientation of the external field. 

III. Results 

By numerically integrating the network equation (Q) 
we computed the macroscopic voltage response function 
of the 3D-network for various different sets of junction 
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Fig. 3. Spherical representation of the macroscopic voltage response 
function of the 3D-network. The time averaged voltage drop Vb 
across the network is plotted in the form V(9, <f>) = Vb eg a as a func- 
tion of the direction eg & of an applied constant external magnetic 
field B ext = S ,ext ^0,0 inducing a flux |/ C xtl= B o,ext a 2 /*0 = 3. 
The system of coordinates on the l.h.s shows the definition of the 
spherical angles 9 and <fi. The origin lies in the center of the cube 
such that Vb = \ V(8, (j>)\. 

and array parameters. In the following we present typ- 
ical results. Fig.[| shows a spherical plot of the voltage 
response function for (3 — 0.5, A = 5, Ig = 4.4 i c and a 
homogeneous magnetic field B ex ^(x, y, z) = B ex ^ eg^ 
which induces a flux |/ ex tl= B Q ex ^ a 2 /$o = 3. Here 
e.g_<p is a vector with unit length pointing in the direc- 
tion given by the spherical coordinates 9 and <j> and the 
origin is put in the center of the 3D-array (c.f. Fig.||). 
The voltage response function is plotted in the form 
V(6, (j>) — Voeg^ = (V x ,V y ,V z ), where Vb is the macro- 
scopic averaged voltage drop across the array measured 
in units of i c r. The eight input and output resistors (c.f. 
Fig.^J) all have the same value so that, corresponding to 
the symmetry of the network, the voltage response func- 
tion shows a periodicity with period tt/2 along each closed 
intersection curve with the x-z-, x-y- and y-x-plane, re- 
spectively, and a reflection symmetry with respect to each 
integer multiple of 7r/4. 

In FigJ|(a) intersection curves in the z-y-plane, i.e., 
= 7r/2, are plotted for 9 = [tt/2, 3tt/2] and different 
values of |/ ex tl- The symmetry of the voltage response 
functions with respect to integer multiples of ir/4 is clearly 
observable and the sensitivity of the network reaches its 
maximum around integer multiples of it/2. For such val- 
ues of 9 the field vector of the external magnetic field lies 
directly perpendicular to one of the faces of the cube and 
the induced currents in general become maximal in the 
direction parallel (and antiparallel) to the bias current. 
In this case the externally induced voltage drop becomes 
also maximal. The number of local maxima and minima 
of the voltage response function is strongly related to the 
applied flux | / Gx t | and can be determined numerically. 

The resolving power of 3D-SQUIDs with respect to 
the strength of external magnetic fields is comparable to 
the resolving power of a conventional 2D-SQUID, i.e. a 
SQUID loop with two junctions. The angular resolution of 
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Fig. 4. (a) Voltage response function V(9, <fi = tt/2) in the y-z- 
plane (x = 0) for 9 = [tt/2, 3-7r/2] and different values of the exter- 
nal applied magnetic field B ex ^ which induces a flux of | / ex t I = 1 
(dotted line), |/ ex tl = 2 (dashed line) and |/ ex tl = 3 (solid line), 
(b) Comnarison of the voltage response function of the 3D-SQUID 
(c.f. Fig.Q) and a conventional 2D-SQUID which consists of a su- 
perconducting loop that contains two junctions. The bias current 
per active junction is for both configurations 1.1 i c . 

3D-SQUIDs, however, is orders of magnitudes better than 
for 2D-SQUIDs. For external magnetic fields inducing a 
flux \fext\ in the order of one, the transfer factor 19^/96*1 
lies around 0.02i c r/rad which is for typical junctions in 
the range of some 10 [iV/rad, such that, e.g., an angle of 
10~ 2 rad should be easily resolvable. The resolving power 
with respect to the strength and the resolving power with 
respect to the direction of an external field can be fur- 
ther increased by using networks which consist of several 
cubes. In this case the resolution powers increase propor- 
tional to the number of cubes. 

Fig.H(b) shows for comparison the angular resolution of 
the 3lXjunction array and a corresponding conventional 
2D-SQUID. Both configurations are considered to have 
identical parameters. If the applied magnetic field lies 
perpendicular to the 2D-SQUID plane, i.e. 6 = tt, the 
device is very insensitive to variations of 9. This insen- 
sitivity, however, has advantages if only the strength of 
the magnetic field is to be measured. In contrast to this, 
the 3D-array's response is very sensitive to variations of 
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the direction of the external magnetic field. Therefore, 
the 3D-network has the advantage that simultaneously 
the strength and the direction of the magnetic field can 
be measured very precisely. 

Three-axis SQUIDs, which consist of three independent 
operating 2D-SQUIDs are not directly comparable with 
3D-SQUIDs. For three-axis SQUIDs the three indepen- 
dent voltage response functions are computationally post- 
processed and the information about strength and orien- 
tation of the magnetic field vector is computationally ex- 
tracted. In contrast to this, the single voltage response 
function of 3D-SQUIDs provides all information at once 
and no further postprocessing is needed. 

For increasing |/ e xtl the number of local maxima and 
minima of the voltage response function grows. Fig]^ 
shows density plots of the response function for a solid 
angle ACl = A9 A(j> with A9 = = 2tt/3, and (a) 
l/extl = ^0 an d (b) l/extl = ^00. The network para- 
meters are the same as for Figs.Q and (^). The minima 
(black dots) and maxima (white areas) form significant 
patterns which are quite similar to optical interference 
patterns. However, according to Eq.(|l|) these patterns are 
a result of nonlinear interactions and a trace of their non- 
linear origin are the complex local structures occurring 
in Fig.||(b). By evaluating the voltage response patterns, 
the strength of the external magnetic field and, up to the 
symmetry redundancy, the orientation of this field can 
be determined very precisely. The symmetry redundancy 
can be removed by an appropriate choice of the input and 
output resistors. In this case, however, the resistors pos- 
sess slightly different values and the voltage response of 
the network becomes much more complicated and more 
difficult to interpret [7]. 



(a) 



(b) 



Fig. 5. Density plot of the voltage response function for a solid 
angle with 9 = [— 7r/3, 7r/3] and <j> = [— 7r/3, 7r/3]. Black dots in- 
dicate minima of the response function and white areas maxima 
(c.f. Fig.^). The external magnetic field induces in (a) a flux of 
|/ ext | = 50andin (b) |/ ext | = 200. 

If the parameters of the array junctions and the di- 
mension of the 3D-network are chosen appropriately, the 
systems are also able to detect time dependent electro- 
magnetic fields with wavelengths which lie typically in 
the mm-range. For time dependent external fields the in- 
duced flux becomes time dependent f cx x = /ext(^) an d 
the currents J r f(t) are induced by the oscillating electric 
field vector of the external wave (c.f. Eq.l). In contrast 
to 2D-configurations, for which the influence of the os- 
cillating magnetic field is in general negligible, for the 
3D-network both field vectors contribute to the response 



function of the array if the lattice spacing a is chosen 
appropriately. Fig.0 shows current-voltage characteris- 
tics of the 3D-network for different fixed directions of the 
external microwave and array parameters = 0.5 and 
A = 5. The microwave is assumed to be linearly polar- 
ized parallel to the bias current and the amplitudes of 
the externally induced flux and of the induced current 
j r £ are |/ e xtl = Jrf = ^he frequency of the exter- 
nal microwave is v = 0.2 i c r in Fig.^(a) and v — 0.8 i c r 
in FigJ|(b), and Ib = Ib/^ is the normalized bias cur- 
rent. In each of both figures the I-V-curves for the mag- 
netic field directions (9, <fi) — (0, 7r/2) (solid curve) and 
(9, (p) — (7r/4, 7r/2) (dashed curve) are plotted, and, for 
comparison, the I-V-curve for negligible influence of the 
magnetic field (dotted curve). 




Fig. 6. Current-voltage characteristics of the 3D-network for a time 
dependent external electromagnetic field. The external microwave is 
linearly polarized parallel to the bias current and the direction of the 
magnetic field vector is in (a) and (b) (0, rf>) = (0, tt/2) (solid curve) 
and (6, <j>) = (tt/A, tt/2) (dashed curve). The frequencies of the 
external microwave are v = 0.2 i c r for (a) and v = 0.8 i c r for (b). 
The amplitudes of the externally induced flux and of the induced 
current are for both figures |/ ex tl = |j'rf| = ^^ ^ or comparison also 
the I-V-curves for negligible influence of the magnetic field (dotted 
curves) are plotted. The dashed and dotted curves in (a) have been 
offset by 0.5 units along the x-axis and in (b) by 0.75 units. The 
inset of (b) shows the splitting of the characteristics for a linearly 
polarized microwave (8, <f>) = (0, tt/2) (solid curve) and a microwave 
that is circularly polarized in the z-x-plane (dashed curve). 

It is clearly observable that the I-V-characteristics de- 



4 



pends on the direction of the incoming wave and that the 
magnetic part is indeed relevant. In general, the distri- 
bution of the Shapiro steps and their width differ signif- 
icantly from those corresponding to a single junction. In 
addition, it can be observed that for some characteristics 
the step plateaus show an unusual fine structure. The 
first, fifth and seventh step on the left characteristics in 
Fig.g(a) show such structures. These fine structures are 
according to Eq.(l) implied by the different contributions 
of the electric and the magnetic field to the network dy- 
namics. We observed that on the Shapiro steps the net- 
work junctions first lock to the electric part of the external 
microwave and then for increasing bias current eventually 
also to the magnetic part. This can change the current 
distributions within the network and implies a slight de- 
crease in the averaged voltage drop and therefore in the 
network frequency. 

The 3D-network dynamics governed by Eq.(l) includes 
all informations about the field variables of the external 
microwave. By a detailed analysis it is therefore possi- 
ble to extract from the macroscopic current and voltage 
response functions of the network the information about 
the phase relationship within the incoming wave. 

By assuming the simple model of a circularly polar- 
ized incoming microwave we can show that the voltage 
response function of the 3D-network is directly affected 
by the hclicity of the external field. The inset of Fig.^(b) 
shows the I-V-characteristics near the critical array cur- 
rent I Cj a for (9, <fi) — (0, 7r/2) for a microwave that is 
linearly polarized in the x-direction (solid curve) and a 
microwave that is circularly polarized in the z-x-plane 
(dashed curve). The splitting between the critical array 
current I c a of the two modes lies in the range of several 
0.01 i c r and the shape of the characteristics differ signifi- 
cantly. This result already indicates that the 3D-network 
can operate phase-sensitively and that the modulations in 
the voltage response function caused by differently polar- 
ized external microwaves should be easily observable in 
experiments. 

IV. Conclusions 

By presenting various theoretical results we have shown 
that 3D-networks of Josephson junctions represent ultra- 
sensitive 3D-SQUIDs which can be useful for a number 
of different applications like magnetic field sensors, video 
detectors or mixers. Especially the capability of the net- 
works to operate direction-sensitively and to detect in- 
coming microwaves phase-sensitively is a novel quality and 
may allow the design of a new generation of superconduct- 
ing devices. 
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